#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
NumPy 实际应用案例
演示NumPy在数据科学、图像处理、信号处理等领域的应用
"""

import numpy as np

def data_preprocessing_example():
    """数据预处理示例"""
    print("=== 数据预处理示例 ===")
    
    # 模拟学生成绩数据
    np.random.seed(42)
    students = 1000
    subjects = 5
    
    # 生成原始成绩数据（包含缺失值）
    scores = np.random.normal(75, 15, (students, subjects))
    scores[scores < 0] = 0
    scores[scores > 100] = 100
    
    # 随机设置一些缺失值
    missing_mask = np.random.random((students, subjects)) < 0.05
    scores[missing_mask] = np.nan
    
    print(f"原始数据形状: {scores.shape}")
    print(f"缺失值数量: {np.sum(np.isnan(scores))}")
    
    # 处理缺失值：用科目平均分填充
    for subject in range(subjects):
        subject_mean = np.nanmean(scores[:, subject])
        scores[np.isnan(scores[:, subject]), subject] = subject_mean
    
    # 数据标准化
    standardized_scores = (scores - np.mean(scores, axis=0)) / np.std(scores, axis=0)
    
    print(f"标准化后均值: {np.mean(standardized_scores, axis=0)}")
    print(f"标准化后标准差: {np.std(standardized_scores, axis=0)}")
    
    # 计算总分和排名
    total_scores = np.sum(scores, axis=1)
    rankings = np.argsort(total_scores)[::-1] + 1  # 降序排名
    
    print(f"最高总分: {np.max(total_scores):.1f}")
    print(f"平均总分: {np.mean(total_scores):.1f}")
    print(f"最低总分: {np.min(total_scores):.1f}")

def image_processing_example():
    """图像处理示例"""
    print("\n=== 图像处理示例 ===")
    
    # 创建模拟图像（灰度）
    height, width = 100, 100
    
    # 创建渐变图像
    x, y = np.meshgrid(np.linspace(0, 1, width), np.linspace(0, 1, height))
    image = np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y)
    image = (image + 1) / 2 * 255  # 归一化到0-255
    
    print(f"图像形状: {image.shape}")
    print(f"像素值范围: {np.min(image):.1f} - {np.max(image):.1f}")
    
    # 添加噪声
    noise = np.random.normal(0, 10, image.shape)
    noisy_image = np.clip(image + noise, 0, 255)
    
    # 简单的均值滤波器（3x3）
    filtered_image = np.zeros_like(image)
    for i in range(1, height-1):
        for j in range(1, width-1):
            filtered_image[i, j] = np.mean(noisy_image[i-1:i+2, j-1:j+2])
    
    # 边缘检测（简单的梯度）
    grad_x = np.gradient(image, axis=1)
    grad_y = np.gradient(image, axis=0)
    gradient_magnitude = np.sqrt(grad_x**2 + grad_y**2)
    
    print(f"梯度幅度范围: {np.min(gradient_magnitude):.1f} - {np.max(gradient_magnitude):.1f}")

def signal_processing_example():
    """信号处理示例"""
    print("\n=== 信号处理示例 ===")
    
    # 生成合成信号
    sampling_rate = 1000
    duration = 2.0
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    
    # 多频率信号
    signal = (np.sin(2 * np.pi * 50 * t) +
              0.5 * np.sin(2 * np.pi * 120 * t) +
              0.3 * np.sin(2 * np.pi * 200 * t))
    
    # 添加噪声
    noise = 0.2 * np.random.randn(len(t))
    noisy_signal = signal + noise
    
    print(f"信号长度: {len(signal)}")
    print(f"采样率: {sampling_rate} Hz")
    print(f"信噪比: {10 * np.log10(np.var(signal) / np.var(noise)):.1f} dB")
    
    # 简单的移动平均滤波
    window_size = 5
    filtered_signal = np.convolve(noisy_signal, 
                                  np.ones(window_size)/window_size, 
                                  mode='same')
    
    # 计算能量
    signal_energy = np.sum(signal**2)
    noise_energy = np.sum(noise**2)
    
    print(f"信号能量: {signal_energy:.1f}")
    print(f"噪声能量: {noise_energy:.1f}")

def financial_analysis_example():
    """金融分析示例"""
    print("\n=== 金融分析示例 ===")
    
    # 模拟股票价格数据
    np.random.seed(42)
    days = 252  # 一年的交易日
    initial_price = 100
    
    # 生成随机收益率（正态分布）
    daily_returns = np.random.normal(0.001, 0.02, days)  # 0.1%均值，2%标准差
    
    # 计算累积价格
    prices = initial_price * np.exp(np.cumsum(daily_returns))
    
    print(f"初始价格: ${initial_price:.2f}")
    print(f"最终价格: ${prices[-1]:.2f}")
    print(f"总收益率: {(prices[-1]/initial_price - 1)*100:.1f}%")
    
    # 计算技术指标
    # 简单移动平均
    window_20 = 20
    sma_20 = np.convolve(prices, np.ones(window_20)/window_20, mode='valid')
    
    # 波动率（20日）
    returns = np.diff(prices) / prices[:-1]
    volatility_20 = np.std(returns[-20:]) * np.sqrt(252) * 100  # 年化波动率
    
    # 最大回撤
    cumulative_max = np.maximum.accumulate(prices)
    drawdown = (prices - cumulative_max) / cumulative_max
    max_drawdown = np.min(drawdown) * 100
    
    print(f"20日移动平均: ${sma_20[-1]:.2f}")
    print(f"年化波动率: {volatility_20:.1f}%")
    print(f"最大回撤: {max_drawdown:.1f}%")

def machine_learning_data_prep():
    """机器学习数据准备示例"""
    print("\n=== 机器学习数据准备示例 ===")
    
    # 生成模拟数据集
    np.random.seed(42)
    n_samples = 1000
    n_features = 5
    
    # 特征矩阵
    X = np.random.randn(n_samples, n_features)
    
    # 目标变量（线性关系 + 噪声）
    true_weights = np.array([1.5, -2.0, 0.5, 3.0, -1.0])
    y = X @ true_weights + 0.1 * np.random.randn(n_samples)
    
    print(f"数据集形状: {X.shape}")
    print(f"目标变量形状: {y.shape}")
    
    # 特征标准化
    X_mean = np.mean(X, axis=0)
    X_std = np.std(X, axis=0)
    X_standardized = (X - X_mean) / X_std
    
    # 数据分割（训练集/测试集）
    split_idx = int(0.8 * n_samples)
    indices = np.random.permutation(n_samples)
    
    train_idx = indices[:split_idx]
    test_idx = indices[split_idx:]
    
    X_train, X_test = X_standardized[train_idx], X_standardized[test_idx]
    y_train, y_test = y[train_idx], y[test_idx]
    
    print(f"训练集大小: {X_train.shape[0]}")
    print(f"测试集大小: {X_test.shape[0]}")
    
    # 简单的线性回归（正规方程）
    # θ = (X^T X)^-1 X^T y
    theta = np.linalg.solve(X_train.T @ X_train, X_train.T @ y_train)
    
    # 预测
    y_pred = X_test @ theta
    
    # 评估
    mse = np.mean((y_test - y_pred)**2)
    r2 = 1 - np.sum((y_test - y_pred)**2) / np.sum((y_test - np.mean(y_test))**2)
    
    print(f"均方误差: {mse:.4f}")
    print(f"R²分数: {r2:.4f}")
    print(f"学习到的权重: {theta}")
    print(f"真实权重: {true_weights}")

def main():
    """运行所有实际应用示例"""
    print("NumPy 实际应用案例演示")
    print("=" * 50)
    
    data_preprocessing_example()
    image_processing_example()
    signal_processing_example()
    financial_analysis_example()
    machine_learning_data_prep()

if __name__ == "__main__":
    main() 